%% fig 3ab orientation selectivity 2 mazes
clear
load(['/media/Big1/15tracks/HeadDirection/head_direction_selectivity_newrand_maze12.mat'],'dsiall2m')



goodind=[dsiall2m(:,1).osi_pval]<0.05 &[dsiall2m(:,2).osi_pval]<0.05;
% goodind=[dsiall2m(:,1).osi_pval]<0.01 ;
% goodind=[dsiall2m(:,1).osi_pval]<0.05 & [dsiall2m(:,1).osi]>0.8&[dsiall2m(:,2).osi_pval]<0.05;
osi1=[dsiall2m(goodind,1).osi]';
dsi1=[dsiall2m(goodind,1).dsi]';
osi2=[dsiall2m(goodind,2).osi]';
dsi2=[dsiall2m(goodind,2).dsi]';
samedir=([dsiall2m(:,1).pref_dir]==[dsiall2m(:,2).pref_dir]);
diffdir=([dsiall2m(:,1).pref_dir]~=[dsiall2m(:,2).pref_dir]);
sameori=([dsiall2m(:,1).pref_ori]==[dsiall2m(:,2).pref_ori]);
diffori=([dsiall2m(:,1).pref_ori]~=[dsiall2m(:,2).pref_ori]);
samedir_rand=([dsiall2m(:,1).pref_dir_rand]==[dsiall2m(:,2).pref_dir_rand]);
diffdir_rand=([dsiall2m(:,1).pref_dir_rand]~=[dsiall2m(:,2).pref_dir_rand]);
sameori_rand=([dsiall2m(:,1).pref_ori_rand]==[dsiall2m(:,2).pref_ori_rand]);
diffori_rand=([dsiall2m(:,1).pref_ori_rand]~=[dsiall2m(:,2).pref_ori_rand]);

ncell=size(dsiall2m,1);
nshuffle=size(samedir_rand,1);

yd=[sum(samedir(goodind)),sum(diffdir(goodind)),sum(samedir(~goodind)),sum(diffdir(~goodind)),...
    sum(samedir_rand(:)),sum(diffdir_rand(:))];
ydp=[sum(samedir(goodind))/sum(goodind),sum(diffdir(goodind))/sum(goodind),...
    sum(samedir(~goodind))/sum(~goodind),sum(diffdir(~goodind))/sum(~goodind),...
    sum(samedir_rand(:))/ncell/nshuffle,sum(diffdir_rand(:))/ncell/nshuffle];

y=[sum(sameori(goodind)),sum(diffori(goodind)),sum(sameori(~goodind)),sum(diffori(~goodind)),...
    sum(sameori_rand(:)),sum(diffori_rand(:))];
yp=[sum(sameori(goodind))/sum(goodind),sum(diffori(goodind))/sum(goodind),...
    sum(sameori(~goodind))/sum(~goodind),sum(diffori(~goodind))/sum(~goodind),...
    sum(sameori_rand(:))/ncell/nshuffle,sum(diffori_rand(:))/ncell/nshuffle];
% f6=figure('position',[328   185   549   503]);
% x=[1,2,4,5,7,8];
% subplot(221);
% bar(x(1:2:end),ydp(1:2:end),0.2);hold on;
% bar(x(2:2:end),ydp(2:2:end),0.2);
% set(gca,'xtick',[1.5,4.5,7.5],'xticklabel',{'sig cell','nonsig cell','shuffle'},'xticklabelrotation',30);
% legend({'Same','Diff'},'location','best');
% title('Direction');
% subplot(222);
% bar(x(1:2:end),yp(1:2:end),0.2);hold on;
% bar(x(2:2:end),yp(2:2:end),0.2);
% set(gca,'xtick',[1.5,4.5,7.5],'xticklabel',{'sig cell','nonsig cell','shuffle'},'xticklabelrotation',30);
% legend({'Same','Diff'},'location','best');
% title('Orientation');
% subplot(223);
% plot_scatter([dsi1,dsi2]);xlabel('DSI(maze1)');ylabel('DSI(maze2)');
% [dsi12_ccr,dsi12_ccp]=corr(dsi1(:),dsi2(:));
% title(sprintf('r=%.2g,p=%.2g',dsi12_ccr,dsi12_ccp));
% subplot(224);
% plot_scatter([osi1,osi2]);xlabel('OSI(maze1)');ylabel('OSI(maze2)');
% [osi12_ccr,osi12_ccp]=corr(osi1(:),osi2(:));
% title(sprintf('r=%.2g,p=%.2g',osi12_ccr,osi12_ccp));
f6=figure('position',[328   185   549   503]);
x=[1,2,4,5,7,8];
subplot(221);
stem(x(1:2:end),ydp(1:2:end));hold on;
stem(x(2:2:end),ydp(2:2:end));
set(gca,'xtick',[1.5,4.5,7.5],'xticklabel',{'sig cell','nonsig cell','shuffle'},'xticklabelrotation',30);
legend({'Same','Diff'},'location','best');
title('Direction');
subplot(222);
stem(x(1:2:end),yp(1:2:end));hold on;
stem(x(2:2:end),yp(2:2:end));
set(gca,'xtick',[1.5,4.5,7.5],'xticklabel',{'sig cell','nonsig cell','shuffle'},'xticklabelrotation',30);
legend({'Same','Diff'},'location','best');
title('Orientation');
subplot(223);
plot_scatter([dsi1,dsi2]);xlabel('DSI(maze1)');ylabel('DSI(maze2)');
[dsi12_ccr,dsi12_ccp]=corr(dsi1(:),dsi2(:));
title(sprintf('r=%.2g,p=%.2g',dsi12_ccr,dsi12_ccp));
subplot(224);
plot_scatter([osi1,osi2]);xlabel('OSI(maze1)');ylabel('OSI(maze2)');
[osi12_ccr,osi12_ccp]=corr(osi1(:),osi2(:));
title(sprintf('r=%.2g,p=%.2g',osi12_ccr,osi12_ccp));

setaxisformalAll(10)

savepdf(f6,fullfile('/media/Big1/15tracks/figures/','OSI_2maze_stem.pdf'));
%% fig 3c s3bc; intersection area. orientation selectivity index -- ori, dir ,location
% 
clear


load(['/media/Big1/15tracks/HeadDirection/Track_intersection4cm/selectivityIndex_ori_dir_loc.mat'],...
    'dir_si_bymean','dir_si_bymean_mean','dir_si_bypeak','dir_si_bypeak_mean','dir_si_out_bymean',...
    'dir_si_out_bymean_mean','dir_si_out_bypeak','dir_si_out_bypeak_mean','loc_si_bymean',...
    'loc_si_bymean_mean','loc_si_bypeak','loc_si_bypeak_mean','si','si_cellid','track_si_bymean',...
    'track_si_bymean_mean','track_si_bypeak','track_si_bypeak_mean','track_si_out_bymean',...
    'track_si_out_bymean_mean','track_si_out_bypeak','track_si_out_bypeak_mean',...
    'ngood0','ngood1','ngood2','pgood2','pgood2m','pyrcell','meanrate','peakrate','npeakrate',...
    'ntrack_1hz','ndir_1hz','ndirtrack_1hz','narea_1hz');
%

f1=figure;
subplot(221)
y=ntrack_1hz;
bar(0:2,y,1,'linestyle','none');set(gca,'ylim',[0,max(y)*1.1]);
xlabel('# of track');ylabel('# of cell~area~dir');
% yyaxis right;bar(0:2,y/sum(y));set(gca,'ylim',[0,max(y/sum(y))*1.1]);
ylabel('prob.');
subplot(222)
y=ndir_1hz;
bar(0:2,y,1,'linestyle','none');set(gca,'ylim',[0,max(y)*1.1]);
xlabel('# of dir');ylabel('# of cell~area~track');
% yyaxis right;bar(0:2,y/sum(y));set(gca,'ylim',[0,max(y/sum(y))*1.1]);
ylabel('prob.');
subplot(223)
y=ndirtrack_1hz;
bar(0:4,y,1,'linestyle','none');set(gca,'ylim',[0,max(y)*1.1]);
xlabel('# of dir/track');ylabel('# of cell~area');
% yyaxis right;bar(0:4,y/sum(y));set(gca,'ylim',[0,max(y/sum(y))*1.1]);
ylabel('prob.');
subplot(224)
y=narea_1hz;
bar(0:6,y,1,'linestyle','none');set(gca,'ylim',[0,max(y)*1.1]);
xlabel('# of area');ylabel('# of cell');
% yyaxis right;bar(0:6,y/sum(y));set(gca,'ylim',[0,max(y/sum(y))*1.1]);
ylabel('prob.');
setaxisformalAll(12,0)

edge=0:.05:1;
f2=figure;
subplot(221);
histogram(track_si_bymean_mean,edge,'Normalization','probability');
title('inside intersection');xlabel('track selectivity index');ylabel('# neuron');
subplot(222);
% histogram(track_si_out_bymean_mean,edge,'Normalization','probability');
% title('outside intersection');xlabel('track selectivity index');ylabel('# neuron');
tmp={track_si_bymean_mean,loc_si_bymean_mean};
[tmpout,grouppair,pval]=test_multigroup_pairwise(cell2mat(tmp));
% plot_errorbar_cell(tmp,{'Orientation','Location'},'median');
% [m1,se1]=mean_se_cell(tmp,'median');
% errorbarformal(1:2,m1,se1,{'Orientation','Location'});
[ydata,gr]=cell2group(tmp');
plot_errorbar_with_data(ydata,gr,{'Orientation','Location'});
sigstar(grouppair,pval);ylabel('Selectivity Index');
set(gca,'xlim',[0.5,2.5],'xticklabelrotation',30)
subplot(223);
histogram(loc_si_bymean_mean,edge,'Normalization','probability');
title('');xlabel('Location selectivity index');ylabel('# neuron');
subplot(224);
tmp={track_si_bymean_mean,track_si_out_bymean_mean,loc_si_bymean_mean};
[tmpout,grouppair,pval]=test_multigroup_pairwise(cell2mat(tmp));
plot_errorbar_cell(tmp,{'IN Track','OUT Track','Loc'},'median');
sigstar(grouppair,pval);ylabel('Selectivity Index');
set(gca,'xlim',[0.5,3.5],'xticklabelrotation',30)
suptitle('by mean');
setaxisformalAll(12)

f3=figure;
subplot(221);
histogram(track_si_bypeak_mean,edge,'Normalization','probability');
title('inside intersection');xlabel('track selectivity index');ylabel('# neuron');
subplot(222);
tmp={track_si_bypeak_mean,loc_si_bypeak_mean};
[tmpout,grouppair,pval]=test_multigroup_pairwise(cell2mat(tmp));
% plot_errorbar_cell(tmp,{'Orientation','Location'},'median');
% [m1,se1]=mean_se_cell(tmp,'median');
% errorbarformal(1:2,m1,se1,{'Orientation','Location'});
[ydata,gr]=cell2group(tmp');
plot_errorbar_with_data(ydata,gr,{'Orientation','Location'});
sigstar(grouppair,pval);ylabel('Selectivity Index');
set(gca,'xlim',[0.5,2.5],'xticklabelrotation',30)
% histogram(track_si_out_bypeak_mean,edge,'Normalization','probability');
% title('outside intersection');xlabel('track selectivity index');ylabel('# neuron');
subplot(223);
histogram(loc_si_bypeak_mean,edge,'Normalization','probability');
title('');xlabel('Location selectivity index');ylabel('# neuron');
subplot(224);
tmp={track_si_bypeak_mean,track_si_out_bypeak_mean,loc_si_bypeak_mean};
[tmpout,grouppair,pval]=test_multigroup_pairwise(cell2mat(tmp));
plot_errorbar_cell(tmp,{'IN Track','OUT Track','Loc'},'median');
sigstar(grouppair,pval);ylabel('Selectivity Index');
set(gca,'xlim',[0.5,3.5],'xticklabelrotation',30)
suptitle('by peak');
setaxisformalAll(12)

%
edge=0:.05:1;
f4=figure('position',[253   104   1022         526]);
subplot(231);
histogram(track_si_bymean,edge,'Normalization','probability');
title('inside intersection');xlabel('track selectivity index');ylabel('# neuron~area');
subplot(232);
histogram(dir_si_bymean,edge,'Normalization','probability');
title('');xlabel('direction selectivity index');ylabel('# neuron~area');
subplot(234);
histogram(track_si_out_bymean,edge,'Normalization','probability');
title('outside intersection');xlabel('track selectivity index');ylabel('# neuron~area');
subplot(235);
histogram(dir_si_out_bymean,edge,'Normalization','probability');
title('');xlabel('direction selectivity index');ylabel('# neuron~area');
subplot(133);
tmp={track_si_bymean,dir_si_bymean,track_si_out_bymean,dir_si_out_bymean};
[tmpout,grouppair,pval]=test_multigroup_rank(tmp);
plot_errorbar_cell(tmp,{'IN Track SI','IN Dir SI','OUT Track SI','OUT Dir SI'});
sigstar(grouppair([1,2,5,6]),pval([1,2,5,6]));
set(gca,'xlim',[0.5,4.5],'xticklabelrotation',30)
setaxisformalAll(12)

suptitle('by mean');

f5=figure('position',[253   104   1022         526]);
subplot(231);
histogram(track_si_bypeak,edge,'Normalization','probability');
title('inside intersection');xlabel('track selectivity index');ylabel('# neuron~area');
subplot(232);
histogram(dir_si_bypeak,edge,'Normalization','probability');
title('');xlabel('direction selectivity index');ylabel('# neuron~area');
subplot(234);
histogram(track_si_out_bypeak,edge,'Normalization','probability');
title('outside intersection');xlabel('track selectivity index');ylabel('# neuron~area');
subplot(235);
histogram(dir_si_out_bypeak,edge,'Normalization','probability');
title('');xlabel('direction selectivity index');ylabel('# neuron~area');
subplot(133);
tmp={track_si_bypeak,dir_si_bypeak,track_si_out_bypeak,dir_si_out_bypeak};
[tmpout,grouppair,pval]=test_multigroup_rank(tmp);
plot_errorbar_cell(tmp,{'IN Track SI','IN Dir SI','OUT Track SI','OUT Dir SI'});
sigstar(grouppair([1,2,5,6]),pval([1,2,5,6]));
set(gca,'xlim',[0.5,4.5],'xticklabelrotation',30)
suptitle('by peak');
setaxisformalAll(12)


savepath='/media/Big1/15tracks/figures/';
savepdf(f1,fullfile(savepath,'intersectionArea_numberOfReponse_tight.pdf'));
savepdf(f2,fullfile(savepath,'intersectionArea_ori_loc_inout_bymean_dot.pdf'));
savepdf(f3,fullfile(savepath,'intersectionArea_ori_loc_inout_bypeak_dot.pdf'));
% savepdf(f4,fullfile(savepath,'intersectionArea_ori_dir_inout_bymean.pdf'));
% savepdf(f5,fullfile(savepath,'intersectionArea_ori_dir_inout_bypeak.pdf'));


%% fig s3de. orientation population vector correlation
clear
load('/media/Big1/15tracks/HeadDirection/Track_intersection4cm/populationVectorCorr.mat','pv'...
    ,'trackpair_str','ccrp_mean_track','ccrp_peak_track','ccrp_mean_dir','ccrp_peak_dir','ccrp_mean_track_shuf'...
    ,'ccrp_peak_track_shuf','ccrp_mean_dir_shuf','ccrp_peak_dir_shuf');

pval_mean=signrank(ccrp_mean_track(:,1),ccrp_mean_dir(:,1));
pval_peak=signrank(ccrp_peak_track(:,1),ccrp_peak_dir(:,1));
f2=figure('position',[ 133   422   857   421]);
subplot(121)
plot_cdf_with_meanse( {ccrp_mean_track(:,1),ccrp_mean_dir(:,1)},{'Orientation','Direction'})
xlabel('PV corr');ylabel('Prop. track pair');title(sprintf('Mean rate. p=%.2g',pval_mean));
subplot(122)
plot_cdf_with_meanse( {ccrp_peak_track(:,1),ccrp_peak_dir(:,1)},{'Orientation','Direction'})
xlabel('PV corr');ylabel('Prop. track pair');title(sprintf('Peak rate. p=%.2g',pval_peak));
setaxisformalAll

pval_mean_ori=ranksum(ccrp_mean_track(:,1),ccrp_mean_track_shuf(:));
pval_mean_dir=ranksum(ccrp_mean_dir(:,1),ccrp_mean_dir_shuf(:));
pval_peak_ori=ranksum(ccrp_peak_track(:,1),ccrp_peak_track_shuf(:));
pval_peak_dir=ranksum(ccrp_peak_dir(:,1),ccrp_peak_dir_shuf(:));

f3=figure('position',[133    64   883   779]);
subplot(221)
plot_cdf_with_meanse( {ccrp_mean_track(:,1),ccrp_mean_track_shuf(:)},{'Data','Shuffle'})
xlabel('PV corr');ylabel('Prop. track pair');title(sprintf('Mean Ori. p=%.2g',pval_mean_ori));
subplot(222)
plot_cdf_with_meanse( {ccrp_peak_track(:,1),ccrp_peak_track_shuf(:,1)},{'Data','Shuffle'})
xlabel('PV corr');ylabel('Prop. track pair');title(sprintf('Peak Ori. p=%.2g',pval_peak_ori));
subplot(223)
plot_cdf_with_meanse( {ccrp_mean_dir(:,1),ccrp_mean_dir_shuf(:)},{'Data','Shuffle'})
xlabel('PV corr');ylabel('Prop. track pair');title(sprintf('Mean Dir. p=%.2g',pval_mean_dir));
subplot(224)
plot_cdf_with_meanse( {ccrp_peak_dir(:,1),ccrp_peak_dir_shuf(:,1)},{'Data','Shuffle'})
xlabel('PV corr');ylabel('Prop. track pair');title(sprintf('Peak Ori. p=%.2g',pval_peak_dir));

setaxisformalAll

%
savepath='/media/Big1/15tracks/figures/';
% savepdf(f2,fullfile(savepath,'pop_vec_ori_dir.pdf'));
% savepdf(f3,fullfile(savepath,'pop_vec_ori_shuffle.pdf'));

%% fig 3d; orientation group pca trajectory 
clear 
load('/media/Big1/15tracks/PowerLaw/ssvar.mat','ssvar')
for irat=5%1:length(ssvar)
    f1=figure('position',[290         105        869         750]);
    ax=[];
    for i=[irat*2-1,irat*2]
        cproj=ssvar(i).cproj;
        runlen=ssvar(i).runlength;
        trackx=ssvar(i).trackx;
        
        load(['/media/Big1/15tracks/',ssvar(i).ratname,'/templtInPlFieldsNEW.mat'],'seqIDx');
        
        horitrack=unique([seqIDx(ismember([seqIDx.HeadDirection],[0,180])).Track]);
        vertitrack=unique([seqIDx(ismember([seqIDx.HeadDirection],[90,270])).Track]);
        
        groupstr={'Horizontal','Vertical'};
        run2plot={horitrack,vertitrack};
        
        %         run2plot={[5,7],[12,14]};
        
        % projection of each track
        xylim=[min(cproj(:,[1,2,3]));max(cproj(:,[1,2,3]))];
        
        nrow=ceil(sqrt(length(runlen)+1));
        xy=[];minlen=min(runlen);
        for irun=1:length(runlen)
            if irun==1
                x=cproj(1:sum(runlen(1:irun)),1,1);
                y=cproj(1:sum(runlen(1:irun)),2,1);
                z=cproj(1:sum(runlen(1:irun)),3,1);
            else
                x=cproj(sum(runlen(1:irun-1))+1:sum(runlen(1:irun)),1,1);
                y=cproj(sum(runlen(1:irun-1))+1:sum(runlen(1:irun)),2,1);
                z=cproj(sum(runlen(1:irun-1))+1:sum(runlen(1:irun)),3,1);
            end
            x=interp1(1:length(x),x,linspace(1,length(x),minlen));
            y=interp1(1:length(y),y,linspace(1,length(y),minlen));
            z=interp1(1:length(z),z,linspace(1,length(z),minlen));
            xy=cat(3,xy,[x(:),y(:),z(:)]);
        end
        mid=floor(minlen/2);
        %     figure;
%         colordata=jet(length(run2plot));
    colordata=[1     0.2     0
     0     1     1];
        ax1=subplot(2,2,(floor(i/irat)-1)*2+1);h1=[];
        for ig=1:length(run2plot)
            irun=ismember(trackx(:,1),run2plot{ig})& trackx(:,2)==1;
            x=squeeze(xy(:,1,irun));y=squeeze(xy(:,2,irun));z=squeeze(xy(:,3,irun));
            h2=plot3(x,y,z,'-','color',colordata(ig,:));hold on; set(h2,'linewidth',2); 
            for il=1:length(h2);h2(il).Color=h2(il).Color*(rand/2+0.5);end
            h1=[h1;h2(1)];
%             h3=plot3(x(end),y(end),z(end),'o','color',colordata(ig,:));
            set(gca,'xlim',xylim(:,1),'ylim',xylim(:,2),'zlim',xylim(:,3),...
                'XGrid','on','YGrid','on','ZGrid','on','xticklabel',[],...
                'yticklabel',[],'zticklabel',[],'TickLength',[0,0]);
        end
        legend([h1],groupstr,'Position',[0.47 0.9 0.09 0.05]);
        xlabel('PC 1');ylabel('PC 2');zlabel('PC 3');
        ax2=subplot(2,2,(floor(i/irat)-1)*2+2);
        for ig=1:length(run2plot)
            irun=ismember(trackx(:,1),run2plot{ig})& trackx(:,2)==1;
            x=squeeze(xy(:,1,irun));y=squeeze(xy(:,2,irun));z=squeeze(xy(:,3,irun));
            
            h2=plot3(x(1,:),y(1,:),z(1,:),'*','color',colordata(ig,:));hold on;
            h3=plot3(x(end,:),y(end,:),z(end,:),'o','color',colordata(ig,:));
            text(x(1,:),y(1,:),z(1,:),num2str(trackx(irun,1)),'fontsize',8);
            %     text(x(mid),y(mid),z(mid),num2str(trackx(irun,1)));
            set(gca,'xlim',xylim(:,1),'ylim',xylim(:,2),'zlim',xylim(:,3),...
                'XGrid','on','YGrid','on','ZGrid','on');
            title([ssvar(i).ratname,'-',num2str(ssvar(i).isShuffle)]);
            %       pause
            %       cla
        end
        ax=[ax;ax1; ax2];
        % figure;
        % plot([squeeze(xy(1,1,:)),squeeze(xy(end,1,:))]',[squeeze(xy(1,2,:)),squeeze(xy(end,2,:))]','o-');
        % text(squeeze(xy(end,1,:)),squeeze(xy(end,2,:)),num2strcell(trackx(:,1)));
    end
    
    Link = linkprop(ax,{'CameraUpVector', 'CameraPosition', 'CameraTarget', 'XLim', 'YLim', 'ZLim'});
    setappdata(gcf, 'StoreTheLink', Link);
    set(gca,'View',[-7.9000 36.4000]);
end
%

savepath='/media/Big1/15tracks/figures/';
savepdf(f1,fullfile(savepath,'ori_population_PCA_trajectory.pdf'));

%% ori by Plfield Corr
clear
load('/media/Big1/15tracks/HeadDirection/head_direction_selectivity_byPlfieldCorr.mat')

edge=0:.02:1;
edgex=edge2x(edge);
n_osi=histcounts(osi,edge,'normalization','probability');hold on;
n_osi_sig=histcounts(osi(osi_pval<0.05),edge);n_osi_sig=n_osi_sig/length(osi);
n_osirand=histcounts(osirand,edge,'normalization','probability');

[h,p_osi]=kstest2(osi,osirand)
addbase=10^-4;
f4=figure('position',[65   585   839   347]);
subplot(121);
scatter(osi+addbase,osi_pval+addbase,'.')
hold on;
plot([addbase,1],0.05*[1,1],'r-');
set(gca,'xscale','linear','yscale','log','box','off','tickdir','out','ylim',...
    [10^-5,2]);
xlabel('Orientation selectivity');
ylabel('p value');
subplot(122);
h1=bar(edgex,n_osi*100,1,'edgecolor','none','facealpha',0.5,'facecolor','b');hold on;
h3=bar(edgex,n_osirand*100,1,'edgecolor','none','facealpha',0.5,'facecolor','r');hold on;
% h2=bar(edgex,n_osi_sig*100,1,'edgecolor','none','facealpha',0.5,'facecolor',[0,0.5,0.5]);
h1=plot(edgex,n_osi*100,'b','linewidth',2);hold on;
% plot(edgex,n_osi_sig*100,':','color',[0,0.5,0.5],'linewidth',2);
h3=plot(edgex,n_osirand*100,'r','linewidth',2);
% set([h1;h2;h3],'linewidth',2)
xlabel('Orientation selectivity');ylabel('%cell');
title(sprintf('%.1f%%significant(%d/%d)p=%.2g',sum(osi_pval<0.05)/length(osi)*100,...
    sum(osi_pval<0.05),length(osi),p_osi));
legend([h1;h3],{'All pyr cell','Shuffle'},'box','off');
box off
setaxisformalAll

savepath='/media/Big1/15tracks/figures/';
savepdf(f4,fullfile(savepath,'ori_byPlfieldCorr.pdf'));
  %% fig s3f. number of subfields
    clear

ratnames={'jsrat1','jsrat2','jsrat3','jsrat4','jsrat5'};

    nf=[];
    
 for irat=1:5
    cd(fullfile('/media/Big1/15tracks/',ratnames{irat}));
% reading and outputting celltype
fid=fopen('celltype','r');
[ct]=textscan(fid,'%s %s');
fclose(fid);
celltype=[strcmp(ct{1},'p'), strcmp(ct{2},'g')];

pyrcell=find(celltype(:,1));


load('RUNTemplateMultiPeak15.mat','templt_seq3')
load templtInPlFieldsNEW.mat seqIDx    
rep1ind=find([seqIDx.Repeat]==1 & [seqIDx.Multitracks]==0);

nf1=zeros(length(pyrcell),length(rep1ind));
for irun=1:length(rep1ind)
    for ic=1:length(pyrcell)
        nf1(ic,irun)=sum(templt_seq3{rep1ind(irun)}(:,1)==pyrcell(ic));
    end
end
nf=cat(1,nf,nf1);
 end
 
 edge=0.5:5.5;edgex=edge2x(edge);
    nfm=histcounts(nf(:),edge);
    nfm=nfm/sum(nfm);
    figure;
    bar(edgex,nfm);
    title(sprintf('mean nfield=%.4g, mean nfield=%.4g\n',mean(nf(:)),...
        mean(nf(nf>0))));
xlabel('# of field'); ylabel('prob.')
disp(sprintf('%.4g %% of place cells have multiple fields',100*mean(sum(nf>1)./sum(nf>0))))
setaxisformal

%   savepdf(gcf,fullfile('/media/Big1/15tracks/figures/','number_of_subfields_distribution.pdf'));
%% fig 3e s3g;remapping of place field.
clear
load('/media/Big1/15tracks/remapfield_difftrack.mat','pfParamall','res');

surround=20; %cm
grid=2;
surroundInd=floor(surround/grid);

pfParam=pfParamall;
%
peakRate_predata=cellfun(@(x) x(surroundInd+1),{pfParam.pf1this})';
peakRate_postdata=cell2mat(cellfun(@(x) x(surroundInd+1,:),{pfParam.pf1other},'UniformOutput',false));
peakRate_postmean=mean(peakRate_postdata,2);
edge=0:.2:20;edgex=edge2x(edge);
peakRate_pre=histcounts(peakRate_predata,edge);
peakRate_post=histcounts(peakRate_postdata,edge);


% newind=cell2mat(cellfun(@(x) max(x,[],1)<minRate, {pfParam.pf1other},'UniformOutput',false))';
% newind=cell2mat(cellfun(@(x) isnan(x), {pfParam.alsoFieldPeak},'UniformOutput',false))';
% oldind=~newind;
newind=cat(1,res.newind);
oldind=cat(1,res.oldind);
difftrackind=cat(1,res.difftrackind);
preind=cat(1,res.preind);
postind=cat(1,res.postind );
paralelind=cat(1,res.paralelind );
orthind=cat(1,res.orthind );
field1ind=cat(1,res.field1ind );
field23ind=cat(1,res.field23ind );
pf_this=cat(1,res.pf_this );
pf_other=cat(1,res.pf_other );

pf_other_pre=cat(1,res.pf_other_pre );
pf_other_post=cat(1,res.pf_other_post);
pf_other_diff=cat(1,res.pf_other_diff );
pf_other_low=cat(1,res.pf_other_low );
pf_other_high=cat(1,res.pf_other_high );
pf_other_same=cat(1,res.pf_other_same );
pf_other_paralel=cat(1,res.pf_other_paralel );
pf_other_orth=cat(1,res.pf_other_orth );
%
pf_other_field1=cat(1,res.pf_other_field1 );
pf_other_field23=cat(1,res.pf_other_field23 );
%

newind_rand=cat(1,res.newind_rand );
oldind_rand=cat(1,res.oldind_rand );
%
pf_other_rand=cat(1,res.pf_other_rand);

pf_other_rand_pre=cat(1,res.pf_other_rand_pre );
pf_other_rand_post=cat(1,res.pf_other_rand_post );
pf_other_rand_diff=cat(1,res.pf_other_rand_diff );
pf_other_rand_low=cat(1,res.pf_other_rand_low );
pf_other_rand_high=cat(1,res.pf_other_rand_high );
pf_other_rand_same=cat(1,res.pf_other_rand_same );
pf_other_rand_paralel=cat(1,res.pf_other_rand_paralel );
pf_other_rand_orth=cat(1,res.pf_other_rand_orth );
%
pf_other_rand_field1=cat(1,res.pf_other_rand_field1 );
pf_other_rand_field23=cat(1,res.pf_other_rand_field23 );
%
%
pf_other_sametrack18=cat(1,res.pf_other_sametrack18 );
pf_other_rand_sametrack18=cat(1,res.pf_other_rand_sametrack18 );

nbin=size(pf_other,2);
data=pf_other_diff-pf_other_rand_diff;
pval_diff=nan(1,nbin); for j=1:nbin;pval_diff(j)=signrank(data(:,j));end
data=pf_other_pre-pf_other_rand_pre;
pval_pre=nan(1,nbin); for j=1:nbin;pval_pre(j)=signrank(data(:,j));end
data=pf_other_post-pf_other_rand_post;
pval_post=nan(1,nbin); for j=1:nbin;pval_post(j)=signrank(data(:,j));end

pval_low=nan(1,nbin); for j=1:nbin;pval_low(j)=ranksum(pf_other_low(:,j),pf_other_rand_low(:,j));end

pval_high=nan(1,nbin); for j=1:nbin;pval_high(j)=ranksum(pf_other_high(:,j),pf_other_rand_high(:,j));end
data=pf_other_same-pf_other_rand_same;
pval_same=nan(1,nbin); for j=1:nbin;pval_same(j)=signrank(data(:,j));end
data=pf_other_paralel-pf_other_rand_paralel;
pval_paralel=nan(1,nbin); for j=1:nbin;pval_paralel(j)=signrank(data(:,j));end
pval_paralel_anova=anova1(data); %-------anova
data=pf_other_orth-pf_other_rand_orth;
pval_orth=nan(1,nbin); for j=1:nbin;pval_orth(j)=signrank(data(:,j));end
pval_orth_anova=anova1(data);
data=pf_other_field1-pf_other_rand_field1;
pval_field1=nan(1,nbin); for j=1:nbin;pval_field1(j)=signrank(data(:,j));end
data=pf_other_field23-pf_other_rand_field23;
pval_field23=nan(1,nbin); for j=1:nbin;pval_field23(j)=signrank(data(:,j));end
data=pf_other_sametrack18-pf_other_rand_sametrack18;
pval_sametrack18=nan(1,nbin); for j=1:nbin;pval_sametrack18(j)=signrank(data(:,j));end

ratname01=cell2mat(cellfun(@(x,y) repmat(x(:),1,length(y)), {pfParam.ratname},{pfParam.otherseqInd},'UniformOutput',false))';
ratname01=mat2cell(ratname01,ones(size(ratname01,1),1),size(ratname01,2));

[grind,gr]=findgroups(ratname01(newind,:));count_new=splitapply(@sum,ones(size(grind)),grind);
[grind,gr]=findgroups(ratname01(oldind,:));count_old=splitapply(@sum,ones(size(grind)),grind);
[grind,gr]=findgroups(ratname01(newind_rand,:));count_new_rand=splitapply(@sum,ones(size(grind)),grind);
[grind,gr]=findgroups(ratname01(oldind_rand,:));count_old_rand=splitapply(@sum,ones(size(grind)),grind);
rigidratio=count_old./(count_new+count_old);
rigidratio_rand=count_old_rand./(count_new_rand+count_old_rand);

f1=figure;
plot_scatter([rigidratio,rigidratio_rand],'o',1,0);
xlabel('proportion of rigid field');ylabel('proportion of rigid field in Shuffle');
setaxisformal

peaks=[pfParam.alsoFieldPeak]; otherpeak=peaks(difftrackind);
peaks_rand=[pfParam.alsoFieldPeakrand]; otherpeak_rand=peaks_rand(difftrackind);

% figure;
% histogram(otherpeak,0:1:21);hold on;
% histogram(otherpeak_rand,0:1:21);

% kruskalwallis
%
% avestr='mean';
% f2=figure('position',[327          80        1063         581]);
% x=[-surroundInd:surroundInd]*grid;
% subplot(341)
% plot_errorbar_v2(pf_this,x,avestr);hold on
% h2=plot_errorbar_v2(pf_other_diff,x,avestr);
% hr=plot_errorbar_v2(pf_other_rand_diff,x,avestr);
% xlabel('Field location');ylabel('Firing rate')
% legend({'Place field','Other track','Shuffle'},'position',[0.01,0.9,0.1,0.075])
% subplot(342)
% h2z=plot_errorbar_v2(pf_other_diff,x,avestr);hold on;
% hrz=plot_errorbar_v2(pf_other_rand_diff,x,avestr);
% sigind=pval_diff<0.05;plot(x(sigind),h2z.YData(sigind),'k*');
% set(h2z,'Color',h2.Color);
% set(hrz,'Color',hr.Color);
% legend([h2z,hrz],{'Other track','Shuffle'},'location','best')
% xlabel('Field location');ylabel('Firing rate')
% subplot(343)
% h2z=plot_errorbar_v2(pf_other_low,x,avestr);hold on;
% hrz=plot_errorbar_v2(pf_other_rand_low,x,avestr);
% sigind=pval_low<0.05;plot(x(sigind),h2z.YData(sigind),'k*');
% set(h2z,'Color',h2.Color);
% set(hrz,'Color',hr.Color);
% xlabel('Field location');ylabel('Firing rate');title('New field');
% subplot(344)
% h2z=plot_errorbar_v2(pf_other_high,x,avestr);hold on;
% hrz=plot_errorbar_v2(pf_other_rand_high,x,avestr);
% sigind=pval_high<0.05;plot(x(sigind),h2z.YData(sigind),'k*');
% set(h2z,'Color',h2.Color);
% set(hrz,'Color',hr.Color);
% xlabel('Field location');ylabel('Firing rate');title('rigid field');
% subplot(345)
% h2z=plot_errorbar_v2(pf_other_same,x,avestr);hold on;
% hrz=plot_errorbar_v2(pf_other_rand_same,x,avestr);
% sigind=pval_same<0.05;plot(x(sigind),h2z.YData(sigind),'k*');
% set(h2z,'Color',h2.Color);
% set(hrz,'Color',hr.Color);
% xlabel('Field location');ylabel('Firing rate');title('Same track');
% subplot(346)
% h2z=plot_errorbar_v2(pf_other_sametrack18,x,avestr);hold on;
% hrz=plot_errorbar_v2(pf_other_rand_sametrack18,x,avestr);
% sigind=pval_sametrack18<0.05;plot(x(sigind),h2z.YData(sigind),'k*');
% set(h2z,'Color',h2.Color);
% set(hrz,'Color',hr.Color);
% xlabel('Field location');ylabel('Firing rate');title('Same track1/track8');
% subplot(347)
% h2z=plot_errorbar_v2(pf_other_pre,x,avestr);hold on;
% hrz=plot_errorbar_v2(pf_other_post,x,avestr);
% h2x=plot_errorbar_v2(pf_other_rand_pre,x,avestr);
% hrx=plot_errorbar_v2(pf_other_rand_post,x,avestr);
% sigind=pval_pre<0.05;plot(x(sigind),h2z.YData(sigind),'k*');
% sigind=pval_post<0.05;plot(x(sigind),hrz.YData(sigind),'k*');
% legend([h2z,hrz,h2x,hrx],{'other track pre','other track post','pre Shuffle','post shuffle'},...
%     'position',[0.57,0.23,0.13,0.1])
% xlabel('Field location');ylabel('Firing rate');title('pre post');
% subplot(348)
% h2z=plot_errorbar_v2(pf_other_paralel,x,avestr);hold on;
% hrz=plot_errorbar_v2(pf_other_orth,x,avestr);
% h2x=plot_errorbar_v2(pf_other_rand_paralel,x,avestr);
% hrx=plot_errorbar_v2(pf_other_rand_orth,x,avestr);
% sigind=pval_paralel<0.05;plot(x(sigind),h2z.YData(sigind),'k*');
% sigind=pval_orth<0.05;plot(x(sigind),hrz.YData(sigind),'k*');
% 
% legend({'parallel other track','orth other track','parallel Shuffle','orth shuffle'},...
%     'position',[0.78,0.23,0.13,0.1])
% xlabel('Field location');ylabel('Firing rate');title('orientation');
% subplot(349)
% h2z=plot_errorbar_v2(pf_other_field1,x,avestr);hold on;
% hrz=plot_errorbar_v2(pf_other_rand_field1,x,avestr);
% set(h2z,'Color',h2.Color);
% set(hrz,'Color',hr.Color);
% sigind=pval_field1<0.05;plot(x(sigind),h2z.YData(sigind),'k*');
% % legend({'other track','shuffle'},'location','best')
% xlabel('Field location');ylabel('Firing rate');title('1st field');
% subplot(3,4,10)
% h2z=plot_errorbar_v2(pf_other_field23,x,avestr);hold on;
% hrz=plot_errorbar_v2(pf_other_rand_field23,x,avestr);
% set(h2z,'Color',h2.Color);
% set(hrz,'Color',hr.Color);
% sigind=pval_field23<0.05;plot(x(sigind),h2z.YData(sigind),'k*');
% xlabel('Field location');ylabel('Firing rate');title('2nd3rd field');
% 
% setaxisformalAll(10);
%=============================lines

avestr='mean';
f3=figure('position',[327          80        1063         581]);
x=[-surroundInd:surroundInd]*grid;
subplot(341)
[m,se]=mean_se(pf_this,avestr);h1=shadedErrorBar(x,m,se,{'r'},0.5);
hold on;
% h2=plot_errorbar_v2(pf_other_diff,x,avestr);
[m,se]=mean_se(pf_other_diff,avestr);h2=shadedErrorBar(x,m,se,{'b'},0.5);
% hr=plot_errorbar_v2(pf_other_rand_diff,x,avestr);
[m,se]=mean_se(pf_other_rand_diff,avestr);hr=shadedErrorBar(x,m,se,{'k'},0.5);
xlabel('Field location');ylabel('Firing rate')
legend([h1.mainLine,h2.mainLine,hr.mainLine],{'Place field','Other track','Shuffle'},'position',[0.01,0.9,0.1,0.075])
subplot(342)
[m,se]=mean_se(pf_other_diff,avestr);h2z=shadedErrorBar(x,m,se,{'b'},0.5);
hold on;
[m,se]=mean_se(pf_other_rand_diff,avestr);hrz=shadedErrorBar(x,m,se,{'k'},0.5);
sigind=pval_diff<0.05;plot(x(sigind),h2z.mainLine.YData(sigind),'k.');
legend([h2z.mainLine,hrz.mainLine],{'Other track','Shuffle'},'location','best')
xlabel('Field location');ylabel('Firing rate')
subplot(343)
[m,se]=mean_se(pf_other_low,avestr);h2z=shadedErrorBar(x,m,se,{'b'},0.5);
hold on;
[m,se]=mean_se(pf_other_rand_low,avestr);hrz=shadedErrorBar(x,m,se,{'k'},0.5);
sigind=pval_low<0.05;plot(x(sigind),h2z.mainLine.YData(sigind),'k.');
xlabel('Field location');ylabel('Firing rate');title('New field');
subplot(344)
[m,se]=mean_se(pf_other_high,avestr);h2z=shadedErrorBar(x,m,se,{'b'},0.5);
hold on;
[m,se]=mean_se(pf_other_rand_high,avestr);hrz=shadedErrorBar(x,m,se,{'k'},0.5);
sigind=pval_high<0.05;plot(x(sigind),h2z.mainLine.YData(sigind),'k.');
xlabel('Field location');ylabel('Firing rate');title('rigid field');
subplot(345)
[m,se]=mean_se(pf_other_same,avestr);h2z=shadedErrorBar(x,m,se,{'b'},0.5);
hold on;
[m,se]=mean_se(pf_other_rand_same,avestr);hrz=shadedErrorBar(x,m,se,{'k'},0.5);
sigind=pval_same<0.05;plot(x(sigind),h2z.mainLine.YData(sigind),'k.');
xlabel('Field location');ylabel('Firing rate');title('Same track');
subplot(346)
[m,se]=mean_se(pf_other_sametrack18,avestr);h2z=shadedErrorBar(x,m,se,{'b'},0.5);
hold on;
[m,se]=mean_se(pf_other_rand_sametrack18,avestr);hrz=shadedErrorBar(x,m,se,{'k'},0.5);
sigind=pval_sametrack18<0.05;plot(x(sigind),h2z.mainLine.YData(sigind),'k.');
xlabel('Field location');ylabel('Firing rate');title('Same track1/track8');
ax1=subplot(347);
[m,se]=mean_se(pf_other_pre,avestr);h2z=shadedErrorBar(x,m,se,{'b'},0.5);hold on;
[m,se]=mean_se(pf_other_rand_pre,avestr);h2zr=shadedErrorBar(x,m,se,{'k'},0.5);
sigind=pval_pre<0.05;plot(x(sigind),h2z.mainLine.YData(sigind),'k.');
xlabel('Field location');ylabel('Firing rate');title('other track pre');
ax2=subplot(348);
[m,se]=mean_se(pf_other_post,avestr);hrz=shadedErrorBar(x,m,se,{'b'},0.5);hold on;
[m,se]=mean_se(pf_other_rand_post,avestr);hrzr=shadedErrorBar(x,m,se,{'k'},0.5);
sigind=pval_post<0.05;plot(x(sigind),hrz.mainLine.YData(sigind),'k.');
xlabel('Field location');ylabel('Firing rate');title('other track post');
ylim1=get(ax1,'ylim');ylim2=get(ax2,'ylim');ylim12=[min([ylim1(1),ylim2(1)]),max([ylim1(2),ylim2(2)])];
set(ax1,'ylim',ylim12);set(ax2,'ylim',ylim12);
ax1=subplot(349);
[m,se]=mean_se(pf_other_paralel,avestr);h2z=shadedErrorBar(x,m,se,{'b'},0.5);hold on;
[m,se]=mean_se(pf_other_rand_paralel,avestr);h2zr=shadedErrorBar(x,m,se,{'k'},0.5);hold on;
sigind=pval_paralel<0.05;plot(x(sigind),h2z.mainLine.YData(sigind),'k.');
xlabel('Field location');ylabel('Firing rate');title('Parallel');
ax2=subplot(3,4,10);
[m,se]=mean_se(pf_other_orth,avestr);hrz=shadedErrorBar(x,m,se,{'b'},0.5);hold on;
[m,se]=mean_se(pf_other_rand_orth,avestr);hrzr=shadedErrorBar(x,m,se,{'k'},0.5);
sigind=pval_orth<0.05;plot(x(sigind),hrz.mainLine.YData(sigind),'k.');
xlabel('Field location');ylabel('Firing rate');title('Orth');
ylim1=get(ax1,'ylim');ylim2=get(ax2,'ylim');ylim12=[min([ylim1(1),ylim2(1)]),max([ylim1(2),ylim2(2)])];
set(ax1,'ylim',ylim12);set(ax2,'ylim',ylim12);
ax1=subplot(3,4,11);
[m,se]=mean_se(pf_other_field1,avestr);h2z=shadedErrorBar(x,m,se,{'b'},0.5);hold on;
[m,se]=mean_se(pf_other_rand_field1,avestr);hrz=shadedErrorBar(x,m,se,{'k'},0.5);
sigind=pval_field1<0.05;plot(x(sigind),h2z.mainLine.YData(sigind),'k.');
xlabel('Field location');ylabel('Firing rate');title('1st field');
ax2=subplot(3,4,12);
[m,se]=mean_se(pf_other_field23,avestr);h2z=shadedErrorBar(x,m,se,{'b'},0.5);hold on;
[m,se]=mean_se(pf_other_rand_field23,avestr);hrz=shadedErrorBar(x,m,se,{'k'},0.5);
sigind=pval_field23<0.05;plot(x(sigind),h2z.mainLine.YData(sigind),'k.');
xlabel('Field location');ylabel('Firing rate');title('2nd3rd field');
 ylim1=get(ax1,'ylim');ylim2=get(ax2,'ylim');ylim12=[min([ylim1(1),ylim2(1)]),max([ylim1(2),ylim2(2)])];
set(ax1,'ylim',ylim12);set(ax2,'ylim',ylim12);
setaxisformalAll(10);

% savepdf(f1,'/media/Big1/15tracks/figures/RigidFieldProportion.pdf')
% % savepdf(f2,'/media/Big1/15tracks/figures/Fieldremapping.pdf')
% savepdf(f3,'/media/Big1/15tracks/figures/FieldremappingLine.pdf')

